Analysis of ligand-receptor interactions was done with the CellPhoneDB package, Efremova et al, 2019. The basic analysis was done in Python as described in the paper, whereafter the data were edited and plotted in R. In case of too many ligand_receptor pairs selection of relevant pairs was made. As in CellPhoneDB some pairs are given as receptor-ligand, this was manually transformed to ligand_receptor by computation in excel. The order of clusters was edited to be the same as in other figures.

Python and CellphoneDB were installed as system resources, please see code.

From the original seurat file a raw count file is generated with ensemble symbol annotations and a metadata file with the cluster annotations. Start is: seu_2021_May.R generated and saved in CMC2.Rmd

#install.packages(c("readxl","writexl")) 

suppressPackageStartupMessages({
library(igraph)
library(ggplot2)
library(biomaRt)
library(Seurat)
library(dplyr)  
library(data.table)
library(reshape2)
library(openxlsx)
library(readxl)
library(xlsx)
library(reshape2)
library(data.table)
#library(reticulate)
library(writexl)
})

path1 <-'/Users/ChristerSylven/Human-Prenatal-Cardiomyocytes/'
path2 <-'/Users/ChristerSylven/'
seu <-  readRDS(paste0(path1, 'seu_2021_May.R'))
#seu
#names(seu[[]])
#table(seu@meta.data$sub_cluster_SAN2) # SAN should be n=8

Idents(object = seu) <- "sub_cluster_SAN2"

seu_Phone <- subset(x = seu, idents = c('M0','M1','M2','M3','M4', 'M5',  'M6',  'M7', 'SAN','0','2','3','4','6','7','9','10','12','14','15','16'))
seu_Phone # all cell clusters except erythocytes

## An object of class Seurat 
## 30646 features across 3371 samples within 2 assays 
## Active assay: SCT (15323 features, 3000 variable features)
##  1 other assay present: RNA
##  2 dimensional reductions calculated: pca, umap

#names(seu_Phone@meta.data)
seu_Phone@meta.data$sub_cluster_SAN2 <- factor(seu_Phone@meta.data$sub_cluster_SAN2)
table(seu_Phone@meta.data$sub_cluster_SAN2)

## 
##  M0  M1  M2  M3  M4  M5  M6  M7 SAN   0   2   3   4   6   7   9  10  12  14  15 
## 190 141 103  79  74  51  57  23   8 652 475 369 340 160 148 122 119 111  76  44 
##  16 
##  29

#count_raw <- seu@raw.data[,data_object@cell.names] 
count_raw <- as.data.frame(GetAssayData(object = seu_Phone, slot = 'counts'))
count_norm <- apply(count_raw, 2, function(x) (x/sum(x))*10000)
counts <- count_norm

ensembl = useMart("ensembl",dataset = "hsapiens_gene_ensembl")
genes  <-  getBM(filters = 'hgnc_symbol',
                 attributes = c('ensembl_gene_id','hgnc_symbol'),
                 values = rownames(counts),
                 mart = ensembl)

counts <- counts[rownames(counts) %in% genes$hgnc_symbol,]

counts <- tibble::rownames_to_column(
  as.data.frame(counts), var = 'hgnc_symbol')

counts <- plyr::join(counts, genes)

counts$hgnc_symbol <- NULL

counts <- cbind(counts[,which(colnames(counts) == 'ensembl_gene_id')], counts)

colnames(counts)[1] <- 'Gene'
counts$ensembl_gene_id <- NULL

metadata <- data.frame(Cell = rownames(seu_Phone@meta.data),
                       cell_type = Idents(object = seu_Phone))

write.table(counts,
            file = './out/counts.txt',
            quote = F,
            col.names = T,
            row.names = F,
            sep = '\t')

write.table(metadata,
            file = './out/metadata.txt',
            quote = F,
            col.names = T,
            row.names = F,
            sep = '\t')

Normalized counts (ensembl_gene_id) and metadata tables are formatted according to Cellphone DB CellphoneDB is then run and heatmap for all ligand_receptor interactions between all clusters generated in the user/out folder:

#library(reticulate)
Sys.setenv(RETICULATE_PYTHON = "/Users/ChristerSylven/Library/Python/3.8/bin/")

# install cellphonedb as a system resource and use path as below
system('python3 -m venv cpdb')
system('source cpdb/bin/activate')
system('pip install cellphonedb')

# to generate means and p values and heatmap between all celltypes use:
# system("cellphonedb method statistical_analysis metadata.txt counts.txt --counts-data=gene_name --iterations=100 --threads=2 --result-precision 2 --output-path=/Users/ChristerSylven/out/")
# system("cellphonedb plot heatmap_plot metadata.txt --pvalues-path=./out/pvalues.txt")
#  alternative:  use cellphone call without path directly in terminal after having installed cellphonedb

Triangular heatmaps of OFT & Conduit CM L-R interactions

To show all ligand_receptor interactions between of non-CM cell types of relevance according to scRNAseq deconvolution on ST maps in OFT (M7) and atrial Conduit (M5) tt containing spots, metafiles for these cell types are saved

#heatmap plot and count_network (table for heatmap) for subgroups
# Metadata for groups
metadata <- read.table('metadata.txt')
colnames(metadata) <- metadata[1,]
metadata <- metadata[-1,]

cc <- c('M7', '2', '7', '6',  '0', '3', '15') # OFT v2
metadata_M7 <- metadata[metadata$cell_type %in% cc,]
metadata_M7[,2]  <- factor(metadata_M7[,2]) # delete cluster w zero cells
write.table(metadata_M7,
            file = './out/metadata_M7.txt',
            quote = F,
           # col.names = T,
            row.names = F,
            sep = '\t')
metadata <- read.table('./out/metadata_M7.txt')
table(metadata_M7[,2])

## 
##   0  15   2   3   6   7  M7 
## 652  44 475 369 160 148  23

system("cellphonedb plot heatmap_plot ./out/metadata_M7.txt --pvalues-path=./out/pvalues.txt")

Following the system command

the OFT heatmap and its count_network data are saved in the out folder

to construct a triangular heatmap of interactions:

Figure 3C

# heatmap with desired order of clusters

count_network <- read.table('./out/count_network.txt', header = TRUE)


#OFT cc <- c('M7', '2', '7', '6',  '0', '3', '15')
count_wide <- setDT(count_network %>% arrange(factor(SOURCE, levels = c('M7', '2', '7', '6',  '0', '3', '15'))))
count_wide <- dcast(count_wide, factor(SOURCE, levels = unique(c('15','3', '0', '6', '7', '2', 'M7'))) ~ factor(TARGET, levels = unique(c('15','3', '0', '6', '7', '2', 'M7'))), value.var = 'count')
count_wide <- count_wide[,-1]
count_wide[lower.tri(-count_wide)] = NA


cc <- c('OFT', 'Epicardium_d',  'Fb_smooth_m', 'Fb_l_vascular', 'Endothelium', 'Fb_skeleton', 'Schwann')

colnames(count_wide) <-  rev(cc)
count_wide$celltype <- rev(cc)
count_wide2 <- count_wide

triangular_plot <- function (count_wide2)
{
count_wide2 <- na.omit(melt(count_wide2, 'celltype', variable='cohort'))
count_wide2$cohort <- factor(count_wide2$cohort)
count_wide2$celltype <- factor(count_wide2$celltype)

#head(count_wide2)
count_wide2$log <- log(count_wide2$value+1)

ggplot(data = count_wide2) +
  aes(x = cohort, y = celltype) +
  ggtitle('') +
  theme_bw() +
  xlab('') +
  ylab('') +
  geom_tile(aes(fill = log), colour = "transparent") +
  scale_fill_gradientn(colours = c("dodgerblue4", 'peachpuff', 'deeppink4')) +
  scale_x_discrete(limits = cc,
                   labels = cc) +
  scale_y_discrete(limits = cc,
                   labels = cc) +
  theme(
        axis.text.x = element_blank(), 
        axis.text.y=element_text(size=12),
        axis.ticks=element_blank(),
        axis.line=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.border=element_blank()) +
annotate("text", x = (1:length(levels(count_wide2$cohort))), y = 1:length(levels(count_wide2$cohort))- 0.75, label = rev(levels(count_wide2$cohort)), size=4, hjust = 1, angle = 90) +
  coord_cartesian(ylim = c(8, -2), clip = "off")
}

triangular_plot(count_wide)

In the same way, create a metadata file for atrial conduit and colocalized cell types and generate count_network and triangular heatmap.

metadata <- read.table('metadata.txt')
colnames(metadata) <- metadata[1,]
metadata <- metadata[-1,]

cc <- c('M5', '0', '2', '6', '7', '12', '15')
metadata_M5<- metadata[metadata$cell_type %in% cc,]
metadata_M5[,2]  <- factor(metadata_M5 [,2]) # delete cluster w zero cells
write.table(metadata_M5,
            file = './out/metadata_M5.txt',
            quote = F,
            #col.names = T,
            row.names = F,
            sep = '\t')

table(metadata_M5[,2])

## 
##   0  12  15   2   6   7  M5 
## 652 111  44 475 160 148  51

system("cellphonedb plot heatmap_plot ./out/metadata_M5.txt --pvalues-path=./out/pvalues.txt")

Figure 4A

count_network <- read.table('./out/count_network.txt', header = TRUE)
count_network[,1] <- as.factor(count_network[,1] )
count_network[,2] <- as.factor(count_network[,2])
count_network$log_count <- log(count_network$count +1 )
count_network$log_count2 <- count_network$log_count

#count_network

#cc <-  c('M5', '6', '2', '3', '12', '4', '7')

# to construct a triangular heatmap of interactions do:
#Conduit + SAN
count_wide <- setDT(count_network %>% arrange(factor(SOURCE, levels = cc)))
count_wide <- dcast(count_wide, factor(SOURCE, levels = unique(c('15','12', '7', '6', '2', '0', 'M5'))) ~ factor(TARGET, levels = unique(c('15','12', '7', '6', '2', '0',  'M5'))), value.var = 'count')
count_wide <- count_wide[,-1]
count_wide[lower.tri(-count_wide)] = NA
#count_wide

cc <-  c('Conduit',  'Endothelium', 'Epicardium_d', 'Fb_l_vascular','Fb_smooth_m', 'Epicardium', 'Schwann')

colnames(count_wide) <- rev(cc)
count_wide$celltype <- rev(cc)
count_wide2 <- count_wide

triangular_plot(count_wide)

Subgroup output

The basic /out/means.txt & /out/pvalues.text are first pruned by deletion of duplicates

Then the relation between integrin and non-integrin ligand_receptor pairs is plotted for interaction between different cardiomyocytes and cardiomyocytes and non-cardiomyocytes.

Integrins / non integrins

S Figure 5A

#subgroup output
means <-  read.delim("/Users/ChristerSylven/out/means.txt", check.names = FALSE)
pvals <- read.delim("/Users/ChristerSylven/out/pvalues.txt", check.names = FALSE)
#pval[pval==0] = 0.0009
#dim(means)
#dim(pvals)

# delete duplicates of L_R / R_L interactoin pairs
pvals0 <- pvals[!duplicated(pvals[,1]),]
#pvals1 <- pvals[unique(pvals[,1]),]

means0 <- means[!duplicated(means[,1]),]
#dim(pvals0)
#dim(means0)
#add L_R / R_L interaction pairs as rowname
rownames(pvals0) <- pvals0[,2]
rownames(means0) <- means0[,2]
#pvals0[1:2,1:10]
#means0[1:2,1:10]
# make dataframe w L_R R_L interaction pairs, secreted, mol 1 rec, mol 2 rec, integrin
pvals_LR <- pvals0[,c(7:9,11)]
pvals_LR <- cbind(rownames(pvals_LR), pvals_LR)
colnames(pvals_LR)[1] <- 'L_R0'

# to find number of interactions and number of integrin interactions for each muscle subtype and subset pvals0 below and get table on number of integrins

is_integrin <- function(cell_interactions) {
  pvals2 <- pvals0[,cell_interactions]
  # keep L_R R_L interaction pair where any value in row has p value < 0.01
  pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
  t2<- sum(pvals_LR[rownames(pvals3),]$is_integrin == 'False')
  t3 <- sum(pvals_LR[rownames(pvals3),]$is_integrin == 'True')
  t1 <- t2 + t3
 # t[1,1] <- t1
  t[1,1] <- t2
  t[1,2] <- t3
  t <<- t # makes dataframe global
}

plot_interactions_integrins <- function(Number) {
  
  ordning <- c("Trabecular", "Compact", "Purkinje", "High G2M", "Exosome-related", 'OFT')
  muscle_muscle <- data.frame(
    Type = rep(c("nonintegrin", "integrin"), each = 6),
    Cell_type = rep((ordning), 2),
    Number = Number)
  
  ggplot(muscle_muscle, aes(x = Cell_type, y = Number)) +
    geom_col(aes(fill = Type), width = 0.7) +
    scale_fill_manual(values = c("#0073C2FF", "#EFC000FF")) + 
    scale_x_discrete(limits = ordning, labels = ordning) +
    ylim(0, 275) +
    theme(
      axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0, size = 32),
          axis.text.y = element_text(size = 32), legend.text =  element_text(size = 32),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
      legend.title = element_text(size=24)
         # axis.text.x=element_blank(),
         # legend.position="none"
          )
}

#theme(axis.line=element_blank(),
     # axis.text.x=element_blank(),
      #axis.text.y=element_blank(),
      #axis.ticks=element_blank(),
      #axis.title.x=element_blank(),
      #axis.title.y=element_blank(),
      
Muscle_Muscle <- data.frame(
  M0 = c( 'M1|M0', 'M2|M0',  'M4|M0',  'M6|M0', 'M7|M0',  'M0|M1', 'M0|M2',  'M0|M4', 'M0|M6', 'M0|M7'),
M1 = c( 'M0|M1', 'M2|M1',  'M4|M1',  'M6|M1', 'M7|M1',  'M1|M0', 'M1|M2',  'M1|M4', 'M1|M6', 'M1|M7'),
M2  = c( 'M0|M2', 'M2|M2',  'M4|M2',  'M6|M2', 'M7|M2',  'M2|M0', 'M2|M1',  'M2|M4', 'M2|M6', 'M2|M7'),
M4 = c( 'M0|M4', 'M1|M4',  'M2|M4',  'M6|M4', 'M7|M4',  'M4|M0', 'M4|M1',  'M4|M2', 'M4|M6', 'M4|M7'),
M6 = c( 'M0|M6', 'M1|M6',  'M2|M6',  'M4|M6', 'M7|M6',  'M6|M0', 'M6|M1',  'M6|M2', 'M6|M4', 'M6|M7'),
M7 =c( 'M0|M7', 'M1|M7',  'M2|M7',  'M4|M7', 'M6|M7',  'M7|M0', 'M7|M1',  'M7|M2', 'M7|M4', 'M7|M6')
)

t <- data.frame()
T <- data.frame()
for (i in 1:6) {
  is_integrin( Muscle_Muscle[,i])
  T <- rbind(T, t)
}


colnames(T) <- c('non_integrins', 'integrins')
T$clusters <- c('M0', 'M1', 'M2', 'M4', 'M6', 'M7')
#T
T2 <- melt(as.data.table(T), id = 'clusters')
#plot_interactions_integrins(T2$value)

ordning <- c("Trabecular", "Compact", "Purkinje", "High G2M", "Exosome-related", 'OFT')
muscle_muscle <- data.frame(
  Type = rep(c("nonintegrin", "integrin"), each = 6),
  Cell_type = rep((ordning), 2),
  Number = T2$value)

p1 <- ggplot(muscle_muscle, aes(x = Cell_type, y = Number)) +
  geom_col(aes(fill = Type), width = 0.7) +
  scale_fill_manual(values = c("#0073C2FF", "#EFC000FF")) + 
  scale_x_discrete(limits = ordning, labels = ordning) +
  ylim(0, 275) +
  theme(
    axis.text.x =  element_blank(),
    axis.text.y = element_text(size = 32), legend.text =  element_text(size = 24),
    axis.title.x=element_blank(),
    axis.title.y=element_blank(),
    legend.title = element_text(size=24),
    # axis.text.x=element_blank(),
    legend.position="none"
  )
 
# Total /Integrin to non CM
# M0 -> M7

Muscle <- data.frame(
  M0 = c('0|M0', '2|M0', '3|M0', '4|M0', '6|M0', '7|M0', '9|M0', '10|M0', '12|M0', '15|M0', '16|M0', 'M0|0', 'M0|2', 'M0|2', 'M0|4', 'M0|6', 'M0|7', 'M0|9', 'M0|10', 'M0|12', 'M0|15', 'M0|16'),
  M1 = c('0|M1', '2|M1', '3|M1', '4|M1', '6|M1', '7|M1', '9|M1', '10|M1', '12|M1', '15|M1', '16|M1', 'M1|0', 'M1|2', 'M1|2', 'M1|4', 'M1|6', 'M1|7', 'M1|9', 'M1|10', 'M1|12', 'M1|15', 'M1|16'),
  M2  = c('0|M2', '2|M2', '3|M2', '4|M2', '6|M2', '7|M2', '9|M2', '10|M2', '12|M2', '15|M2', '16|M2', 'M2|0', 'M2|2', 'M2|2', 'M2|4', 'M2|6', 'M2|7', 'M2|9', 'M2|10', 'M2|12', 'M2|15', 'M2|16'),
  M4 = c('0|M4', '2|M4', '3|M4', '4|M4', '6|M4', '7|M4', '9|M4', '10|M4', '12|M4', '15|M4', '16|M4', 'M4|0', 'M4|2', 'M4|2', 'M4|4', 'M4|6', 'M4|7', 'M4|9', 'M4|10', 'M4|12', 'M4|15', 'M4|16'),
  M6 = c('0|M6', '2|M6', '3|M6', '4|M6', '6|M6', '7|M6', '9|M6', '10|M6', '12|M6', '15|M6', '16|M6', 'M6|0', 'M6|2', 'M6|2', 'M6|4', 'M6|6', 'M6|7', 'M6|9', 'M6|10', 'M6|12', 'M6|15', 'M6|16'),
  M7 =c('0|M7', '2|M7', '3|M7', '4|M7', '6|M7', '7|M7', '9|M7', '10|M7', '12|M7', '15|M7', '16|M7', 'M7|0', 'M7|2', 'M7|2', 'M7|4', 'M7|6', 'M7|7', 'M7|9', 'M7|10', 'M7|12', 'M7|15', 'M7|16')
)

t <- data.frame()
T <- data.frame()
for (i in 1:6) {
  is_integrin( Muscle[,i])
  T <- rbind(T, t)
}
colnames(T) <- c('non-integrins', 'integrins')
T$clusters <- c('M0', 'M1', 'M2', 'M4', 'M6', 'M7')
#T
T2 <- melt(as.data.table(T), id = 'clusters')
#plot_interactions_integrins(T2$value)

muscle_muscle <- data.frame(
  Type = rep(c("nonintegrin", "integrin"), each = 6),
  Cell_type = rep((ordning), 2),
  Number = T2$value)

p2 <- ggplot(muscle_muscle, aes(x = Cell_type, y = Number)) +
  geom_col(aes(fill = Type), width = 0.7) +
  scale_fill_manual(values = c("#0073C2FF", "#EFC000FF")) + 
  scale_x_discrete(limits = ordning, labels = ordning) +
  ylim(0, 275) +
  theme(
    #axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0, size = 32),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0, size = 24),
    axis.text.y = element_text(size = 24), legend.text =  element_text(size = 24),
    axis.title.x=element_blank(),
    axis.title.y=element_blank(),
    legend.title = element_text(size=24),
    # axis.text.x=element_blank(),
    legend.position="bottom"
  )

p1/p2

To find Ligand_receptor dotplot according to cellphone DB and then edit results for

selected top LR pairs
rearrange so that pairs are shown as L-R and not R_L and make dotplot.

Plot as out/cellphoneDB plot and decide which pairs are top pairs
If out/cellphoneDB only show limited interactions as for M6, M7 analysis and ordering can be made directly manually.

OFT and its interactions

system(‘cellphonedb plot dot_plot –rows in/OFT_pvals2_Nov_rows.txt –columns in/OFT_Nov.txt –output-name OFT_Nov.png’)

OFT_Nov.png saved in out folder with only 27 interactions, no selection is needed

#,2, 7, 6,  0, 3, 15 
OFT <- c('2|M7', '7|M7', '6|M7', '0|M7', '3|M7', '15|M7', 'M7|2', 'M7|7', 'M7|6', 'M7|0', 'M7|3', 'M7|15')
write.table(OFT, file = '/Users/ChristerSylven/in/OFT_Nov.txt', sep = '\t', quote = F, row.names = F, col.names = F)

pvals2 <- pvals0[,OFT]
# to minimze number of interactions, filtered to those with p<0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
means3 <- means0[rownames(pvals3),OFT]
dim(pvals3) # 27 12

## [1] 27 12

write.table(rownames(pvals3), '/Users/ChristerSylven/in/OFT_pvals2_Nov_rows.txt', sep = '\t', quote = F, row.names = F, col.names = F)
OFT <- c('2|M7', '7|M7', '6|M7', '0|M7', '3|M7', '15|M7', 'M7|2', 'M7|7', 'M7|6', 'M7|0', 'M7|3', 'M7|15')


pvals4_calc <- function(name, LRmrna, LRcells )
{
  pvals2 <- pvals0[,LRcells] #pval0 duplicates eliminated above
  pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
  pvals3 <- pvals3[LRmrna,]
  means3 <- means0[rownames(pvals3), LRcells]
  pvals_LR_rows <- pvals_LR[rownames(pvals3),]
  pvals4 <- cbind(pvals3,pvals_LR_rows)
  means4 <- cbind(means3,pvals_LR_rows)
  writexl::write_xlsx(pvals4, paste0("pvals4_d_", name, ".xlsx"))
 writexl::write_xlsx(means4, paste0("means4_d_", name, ".xlsx"))
}

pvals4_calc('OFT', rownames(pvals3), OFT) # 27 pairs of which 12 are integrins

Edit in excel. Change RL to LR into L_R column. Change values accordingly.

For those changed correct mean molecule1/molecule2 in excel to molecule2/molecule1 (1/value)

Save as b for both corrected pvals and means files with all pairs as L_R.

Figure 3D

pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_OFT_b.xlsx'))
means4b <- read_xlsx(paste0(path1, 'means4_d_OFT_b.xlsx'))
# edit which pairs and order
#namn_OFT <- c(8:10,12,13:23,27, 24:25,28:30)

# function for manual dotplot, order and edited L_R pairs given by namn, default all LR pairs in 4b files
dotplot_LR <- function(pvals4b, means4b, L_Rs, length) # length <- length(OFT) etc
{
L_R  <- length + 6
pvals4b <- as.data.frame(pvals4b,  row.names=NULL)
means4b <- as.data.frame(means4b)
rownames(pvals4b) <- pvals4b[,L_R]
pvals4b <- pvals4b[,c(L_R, 1:length)]
pvals4b[pvals4b==0] = 0.0009

means4b <- means4b[,c(L_R, 1:length)]
means4b[means4b==0] = 1

pvals4b$L_R <- factor(pvals4b$L_R, levels = pvals4b$L_R )
means4b$L_R <- factor(means4b$L_R, levels = means4b$L_R )

pvals5 <- reshape2::melt(pvals4b, id=c('L_R'))
means5 <- reshape2::melt(means4b, id = c('L_R'))

level1 <- pvals5[1:dim(pvals4b)[1],1]
plot.data = cbind(pvals5,log2(means5[,3]))

colnames(plot.data) = c('pair', 'clusters', 'pvalue', 'mean')

#my_palette <- colorRampPalette(c("black", "blue", 'grey',"yellow", "red"), alpha=TRUE)(n=399)
my_palette <- colorRampPalette(c("black", "yellow","red"), alpha=TRUE)(n=399)

p <- ggplot2::ggplot(plot.data,aes(x=clusters,y=pair)) +
  geom_point(aes(size=-log10(pvalue),color=mean)) +
  scale_color_gradientn('Log2 mean (Ligand_receptor)', colors=my_palette) +
  #scale_y_discrete(limits = level1, labels = level1)+
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        axis.text=element_text(size=12, colour = "black"),
        axis.text.x = element_text(angle = 90, face = 'italic', hjust = 1, vjust=1),
        axis.text.y = element_text(size=12,  colour = "black"),
        axis.title=element_blank(),
        panel.border = element_rect(size = 0.7, linetype = "solid", colour = "black"),
        legend.position="top",
        plot.margin = unit(c(0.1,1,0,2), "cm")) +
        scale_x_discrete(
          labels = L_Rs) +
coord_flip()

#ggsave("/Users/ChristerSylven/Desktop/CMC_ms/GitHub_code/non_Cond_SAN_CM_L_R.pdf", p, height = 10, width = 12,  device = "pdf")

p
}

L_Rs_OFT <- c("Epicardium_d~OFT",
"Fb_smooth_m~OFT",
"Fb_l_vascular~OFT",
"Endothelium~OFT",
"Fb_skeleton~OFT",
"Schwann~OFT", 
"OFT~Epicardium_d",
"OFT~Fb_smooth_m",
"OFT~Fb_l_vascular",
"OFT~Endothelium",
"OFT~Fb_skeleton",
"OFT~Schwann")

dotplot_LR(pvals4b, means4b, L_Rs_OFT, length(OFT))

Schwann cell interactions

Figure S5B

non-integrin L-R interactions between Schwann progenitor cells and non CM celltypes colocalized in the outflow tract

Schwann 15 Epicardium_d 2, Fb_smooth-m 10, Fb_l_vascular 6, Endothelium 9, Fb_skeleton 3 All non-integrins interactions are shown.

Schwann <- c('2|15', '10|15', '6|15', '9|15', '3|15',   '15|2', '15|10', '15|6', '15|9', '15|3')

pvals2 <- pvals0[,Schwann] #153

# only non-integrins
pvals0_1 <- pvals0[which( pvals0$is_integrin == 'False'),]
pvals2 <- pvals0_1[,Schwann]

# to minimze number of interactions, filtered to those with p<0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
means3 <- means0[rownames(pvals3),Schwann]
dim(pvals3) # 61 10

## [1] 61 10

#add characteristics of each L-R
pvals_LR_Schwann <-pvals_LR[rownames(pvals3),]

  pvals4 <- cbind(pvals3,pvals_LR_Schwann)
  means4 <- cbind(means3,pvals_LR_Schwann)

#Schwann_non_integrin_rows2 L-R pairs ordered manually in Excel intocategories

Schwann_non_integrin_rows2 <- read.delim(paste0(path1,'data/Schwann_non_integrin_rows2.txt'), header = F)

#L-R pairs in Excel  manually reordered into categories  and saved as xlxx file
pvals4_reorder <- pvals4[match(Schwann_non_integrin_rows2[,1], rownames(pvals4)),]
means4_reorder <- means4[match(Schwann_non_integrin_rows2[,1], rownames(means4)),]

writexl::write_xlsx(pvals4_reorder, paste0(path1,"/pvals4_", 'Schwann', ".xlsx"))
writexl::write_xlsx(means4_reorder, paste0(path1,"/means4_", 'Schwann', ".xlsx"))

# edit in excel R_L to L_R, add column L_R with correct L_R names, invert transformed values in means, save as:pvals4_Schwann_d.xlsx, means4_Schwann_d.xlsx

L_Rs_Schwann <- c(
"Epicardium_d~Schwann",
"Fb_smooth_m~Schwann",
"Fb_l_vascular~Schwann",
"Endothelium~Schwann",
"Fb_skeleton~Schwann",
"Schwann~Epicardium_d", 
"Schwann~Fb_smooth_m",
"Schwann~Fb_l_vascular",
"Schwann~Endothelium",
"Schwann~Fb_skeleton"
)

pvals4b <- read_xlsx(paste0(path1,'pvals4_Schwann_d.xlsx'))
means4b <- read_xlsx(paste0(path1,'means4_Schwann_d.xlsx'))


dotplot_LR(pvals4b, means4b, L_Rs_Schwann, length(Schwann))

Exosome-enriched interactions

EXO <- c('2|M6', '7|M6', '6|M6', '0|M6', '3|M6', '15|M6', 'M6|2', 'M6|7', 'M6|6', 'M6|0', 'M6|3', 'M6|15')
write.table(EXO, file  = './in/EXO_Nov.txt', sep = '\t', quote = F, row.names = F, col.names = F)
pvals2 <- pvals0[,EXO]
# to minimze number of interactions, filtered to those with p<0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
means3 <- means0[rownames(pvals3),EXO]
dim(pvals3) # 31 12

## [1] 32 12

write.table(rownames(pvals3), './in/EXO_pvals2_Nov_rows.txt', sep = '\t', quote = F, row.names = F, col.names = F)


pvals4_calc('EXO', rownames(pvals3), EXO) # 32 pairs of which 11 are integrins

system(‘cellphonedb plot dot_plot –rows in/EXO_pvals2_Nov_rows.txt –columns in/EXO_Nov.txt –output-name EXO_Nov.png’)

S Figure 5C

# edit in excel. Change RL to LR into L_R column. Change values accordingly. Save as b for both files.
# Correct mean molecule1/molecule2 in excel with  to molecule2/molecule1 (1/value)
# 
#corrected files with all pairs as L_R saved as b
pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_EXO_b.xlsx'))
means4b <- read_xlsx(paste0(path1, 'means4_d_EXO_b.xlsx'))

L_Rs_EXO <- c("Epicardium_d~EXO",
"Fb_smooth_m~EXO",
"Fb_l_vascular~EXO",
"Endothelium~EXO",
"Fb_skeleton~EXO",
"Schwann~EXO", 
"EXO~Epicardium_d",
"EXO~Fb_smooth_m",
"EXO~Fb_l_vascular",
"EXO~Endothelium",
"EXO~Fb_skeleton",
"EXO~Schwann")
# Edited L_R order:
# namn_exo <- c(2:19,23:25,22, 28:29,20:21,27,30,32:33,26,31)
dotplot_LR(pvals4b, means4b, L_Rs_EXO, length(EXO))

Atrial conduit interactions

# calculate number of L_Rs and integrin L_Rs in Conduit CMs
# 
Conduit_lig<-  c('M5|0', 'M5|2', 'M5|6', 'M5|7', 'M5|12', 'M5|15')
pvals2 <- pvals0[,Conduit_lig]
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3) # 102

LR_cond <- dim(pvals3)[1]

pvals0_1 <- pvals0[which( pvals0$is_integrin == 'False'),]
pvals2_1 <- pvals0_1[,Conduit_lig]
pvals3_1 <- pvals2_1[apply(pvals2_1, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3_1) # 70
LR_cond_non_integrin <- dim(pvals3_1)[1]

Conduit_rec<-  c('0|M5', '2|M5', '6|M5', '7|M5', '12|M5', '15|M5')
pvals2 <- pvals0[,Conduit_rec]
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3) # 124

LR_cond <- LR_cond + dim(pvals3)[1]

pvals0_1 <- pvals0[which( pvals0$is_integrin == 'False'),]
pvals2_1 <- pvals0_1[,Conduit_rec]
pvals3_1 <- pvals2_1[apply(pvals2_1, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3_1)

LR_cond_non_integrin <- LR_cond_non_integrin + dim(pvals3_1)[1]

# VCAM1 Ig adhesion molecule viewed to have non integrin function and added to non integrins

cat('Number of Conduit CM total L-R interactions: ', LR_cond, ' and integrin interactions: ', (LR_cond - LR_cond_non_integrin - 2),
'\nVCAM1 Ig adhesion molecule viewed to have non integrin function and added to non integrins.',
'\nFigure shows L-R interactions for non integrins and for non collagen integrins.')

## Number of Conduit CM total L-R interactions:  226  and integrin interactions:  72 
## VCAM1 Ig adhesion molecule viewed to have non integrin function and added to non integrins. 
## Figure shows L-R interactions for non integrins and for non collagen integrins.

# edited order of Conduit L-R and R-l receptor pairs. Non collagen integrins added (n=10)
COND_lig_rec <- read_xlsx(paste0(path1, 'data/COND_lig_rec.xlsx')) # Edited order of non collagen L_R pairs for Counduit (L-R) and Conduit (R_L)

CONDUIT_lig <- as.data.frame(na.omit(COND_lig_rec[,1]))
pvals4_calc('Cond_lig2', CONDUIT_lig[,1], Conduit_lig)

CONDUIT_rec <- as.data.frame(na.omit(COND_lig_rec[,2]))
pvals4_calc('Cond_rec2', CONDUIT_rec[,1], Conduit_rec)

# edit in excel. Change RL to LR into L_R column. Change values accordingly. Save as b for both pvalue4 and mean4  files.
# Correct mean molecule1/molecule2 in excel with  to molecule2/molecule1 (1/value)
# 
#corrected files with all pairs as L_R saved as b 
# 

L_Rs_nonCM_COND_lig <- c(
"Conduit~Endothelium",
"Conduit~Epicardium_d",
"Conduit~Fb_l_vascular",
"Conduit~Fb_smooth_m",
"Conduit~Epicardium",
"Conduit-Schwann")

L_Rs_nonCM_COND_rec <- c(
"Endothelium~Conduit",
"Epicardium_d~Conduit",
"Fb_l_vascular~Conduit",
"Fb_smooth_m~Conduit",
"Epicardium~Conduit",
"Schwann~Conduit")

Figure 4 B & C

pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_COND_lig_b.xlsx'))#87
means4b <- read_xlsx(paste0(path1, 'means4_d_COND_lig_b.xlsx'))
d1 <- dotplot_LR(pvals4b, means4b,L_Rs_nonCM_COND_lig, length(Conduit_lig))


pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_COND_rec_b.xlsx')) #75 87+75 = 162  CONDUIT_lig (66) + CONDUIT_rec (80) = 146 + 16? 11 non coll integrins
means4b <- read_xlsx(paste0(path1, 'means4_d_COND_rec_b.xlsx'))

d2 <- dotplot_LR(pvals4b, means4b,L_Rs_nonCM_COND_rec, length(Conduit_rec))

d1 / d2

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] writexl_1.4.0      xlsx_0.6.5         readxl_1.4.0       openxlsx_4.2.5    
##  [5] reshape2_1.4.4     data.table_1.14.2  dplyr_1.0.9        sp_1.5-0          
##  [9] SeuratObject_4.1.0 Seurat_4.1.1       biomaRt_2.52.0     ggplot2_3.3.6     
## [13] igraph_1.3.4      
## 
## loaded via a namespace (and not attached):
##   [1] BiocFileCache_2.4.0    plyr_1.8.7             lazyeval_0.2.2        
##   [4] splines_4.2.0          listenv_0.8.0          scattermore_0.8       
##   [7] GenomeInfoDb_1.32.2    digest_0.6.29          htmltools_0.5.3       
##  [10] fansi_1.0.3            magrittr_2.0.3         memoise_2.0.1         
##  [13] tensor_1.5             cluster_2.1.3          ROCR_1.0-11           
##  [16] globals_0.15.1         Biostrings_2.64.0      matrixStats_0.62.0    
##  [19] spatstat.sparse_2.1-1  prettyunits_1.1.1      colorspace_2.0-3      
##  [22] blob_1.2.3             rappdirs_0.3.3         ggrepel_0.9.1         
##  [25] xfun_0.31              crayon_1.5.1           RCurl_1.98-1.7        
##  [28] jsonlite_1.8.0         progressr_0.10.1       spatstat.data_2.2-0   
##  [31] survival_3.3-1         zoo_1.8-10             glue_1.6.2            
##  [34] polyclip_1.10-0        gtable_0.3.0           zlibbioc_1.42.0       
##  [37] XVector_0.36.0         leiden_0.4.2           future.apply_1.9.0    
##  [40] BiocGenerics_0.42.0    abind_1.4-5            scales_1.2.0          
##  [43] DBI_1.1.3              spatstat.random_2.2-0  miniUI_0.1.1.1        
##  [46] Rcpp_1.0.9             viridisLite_0.4.0      xtable_1.8-4          
##  [49] progress_1.2.2         reticulate_1.25-9000   spatstat.core_2.4-4   
##  [52] bit_4.0.4              stats4_4.2.0           htmlwidgets_1.5.4     
##  [55] httr_1.4.3             RColorBrewer_1.1-3     ellipsis_0.3.2        
##  [58] ica_1.0-3              farver_2.1.1           rJava_1.0-6           
##  [61] pkgconfig_2.0.3        XML_3.99-0.10          uwot_0.1.11           
##  [64] deldir_1.0-6           sass_0.4.2             dbplyr_2.2.1          
##  [67] utf8_1.2.2             labeling_0.4.2         tidyselect_1.1.2      
##  [70] rlang_1.0.4            later_1.3.0            AnnotationDbi_1.58.0  
##  [73] cellranger_1.1.0       munsell_0.5.0          tools_4.2.0           
##  [76] cachem_1.0.6           cli_3.3.0              generics_0.1.3        
##  [79] RSQLite_2.2.15         ggridges_0.5.3         evaluate_0.15         
##  [82] stringr_1.4.0          fastmap_1.1.0          goftest_1.2-3         
##  [85] yaml_2.3.5             knitr_1.39             bit64_4.0.5           
##  [88] fitdistrplus_1.1-8     zip_2.2.0              purrr_0.3.4           
##  [91] RANN_2.6.1             KEGGREST_1.36.3        nlme_3.1-158          
##  [94] pbapply_1.5-0          future_1.27.0          mime_0.12             
##  [97] xml2_1.3.3             compiler_4.2.0         rstudioapi_0.13       
## [100] plotly_4.10.0          filelock_1.0.2         curl_4.3.2            
## [103] png_0.1-7              spatstat.utils_2.3-1   tibble_3.1.8          
## [106] bslib_0.4.0            stringi_1.7.8          highr_0.9             
## [109] rgeos_0.5-9            lattice_0.20-45        Matrix_1.4-1          
## [112] vctrs_0.4.1            pillar_1.8.0           lifecycle_1.0.1       
## [115] spatstat.geom_2.4-0    lmtest_0.9-40          jquerylib_0.1.4       
## [118] RcppAnnoy_0.0.19       cowplot_1.1.1          bitops_1.0-7          
## [121] irlba_2.3.5            httpuv_1.6.5           patchwork_1.1.1       
## [124] R6_2.5.1               promises_1.2.0.1       KernSmooth_2.23-20    
## [127] gridExtra_2.3          IRanges_2.30.0         parallelly_1.32.1     
## [130] codetools_0.2-18       MASS_7.3-58            assertthat_0.2.1      
## [133] xlsxjars_0.6.1         withr_2.5.0            sctransform_0.3.3     
## [136] S4Vectors_0.34.0       GenomeInfoDbData_1.2.8 mgcv_1.8-40           
## [139] parallel_4.2.0         hms_1.1.1              rpart_4.1.16          
## [142] grid_4.2.0             tidyr_1.2.0            rmarkdown_2.14        
## [145] Rtsne_0.16             Biobase_2.56.0         shiny_1.7.2

CMC Cellphone analysis

Christer Sylvén

22-11-20